AEP.error.input_msg <- "The correct syntax is 'AEP(frame = f, value=v)' where 'f' is of type 'aep.frame' and 'v' is numeric, \n
  or 'AEP(data = d, alpha = a, nb_iteration = n, value=v)' where 'data' is a 'data.frame' and 'alpha', 'nb_iteration' and 'value' are numeric."

# Define an AEP class
setClass("AEP.frame", representation(
	data = "data.frame",
	pobs_data = "matrix",
	nb_obs = "numeric",
	nb_var = "numeric",
	EmpiricalCDF = "list",
	ECDF = "function",
	N = "numeric",	
	binary_matrix = "array",
	sum_ones = "array",
	Hmultiplier = "array",
	alpha = "numeric",
	nb_iteration = "numeric",
	m = "array"
	)
)

setMethod("initialize", "AEP.frame", function(.Object, data, alpha, nb_iteration) {
  if(nargs() < 3) stop("missing arguments")     
	if (is.data.frame(data) == FALSE) stop("the first parameter must be of type 'data.frame'") 			
	.Object@data <- data
	.Object@nb_obs <- dim(data)[1]
  .Object@nb_var <- dim(data)[2]  
	.Object@alpha <- alpha
	if (alpha < 1 / .Object@nb_var || alpha >= 1) stop("'alpha' must be in [1/d, 1)") 			
	.Object@nb_iteration <- nb_iteration
	.Object@N <- 2^.Object@nb_var
	.Object@binary_matrix <- array(0, c(.Object@N, .Object@nb_var))			
	.Object@sum_ones <- array(0, c(.Object@N))
	.Object@Hmultiplier <- array(0, c(.Object@N))
			
	EmpiricalCDF <- NULL
	#Compute the ECDF of each variable
	for (i in 1:.Object@nb_var) {  
		.Object@EmpiricalCDF <- c(.Object@EmpiricalCDF, ecdf(.Object@data[,i]))
	}
	.Object@ECDF <- function(i, x) .Object@EmpiricalCDF[[i]](x)

  # Compute the pobs of data
  .Object@pobs_data <- pobs(.Object@data)
  #.Object@pobs_data <- array(0, c(.Object@nb_obs, .Object@nb_var))
  #for (i in 1:.Object@nb_obs) {
  #  for (j in 1:.Object@nb_var) {
  #    .Object@pobs_data[i,j] <- .Object@ECDF(j, .Object@data[i,j]) 
  #  }      
  #}
  #.Object@pobs_data <- as.matrix(.Object@pobs_data)
  
  #Initialize the matrix of binary and the sum of ones in each line
	for (i in 2:.Object@N) {
		.Object@binary_matrix[i,] <- .Object@binary_matrix[i-1,]
		j <- 1
		.Object@binary_matrix[i, j] <- .Object@binary_matrix[i, j] + 1
		while (.Object@binary_matrix[i, j] > 1) {
			.Object@binary_matrix[i, j] <- 0
			j <- j + 1
			.Object@binary_matrix[i, j] <- .Object@binary_matrix[i, j] + 1
		}
		.Object@sum_ones[i] <- sum(.Object@binary_matrix[i,])
	}

	#Initialize the multiplier of H()
	for (i in 1:.Object@N) {
		if ((.Object@nb_var  -  .Object@sum_ones[i]) %% 2 == 0) .Object@Hmultiplier[i] <- 1
		else  .Object@Hmultiplier[i] <- -1
	}

	# Initialize the vector m
	.Object@m <- array(0, c(.Object@N))
	.Object@m[1] <- 0
	for (j in 2:.Object@N) {
		if (.Object@sum_ones[j] < 1 / .Object@alpha) {
			if((1 + .Object@sum_ones[j]) %% 2 == 0) .Object@m[j] <- 1
			else .Object@m[j] <- -1
		}
		else if (.Object@sum_ones[j] == 1 / .Object@alpha) .Object@m[j] <- 0
		else {
			if((.Object@nb_var + 1 - .Object@sum_ones[j]) %% 2 == 0) .Object@m[j] <- 1
			else .Object@m[j] <- -1
		}
	}
	.Object
})

setMethod(f = "show", signature = "AEP.frame", definition = function(object) {
  cat("Nb observations =", object@nb_obs, "\nNb variables =", object@nb_var, "\n")
	print(object@data)
})

# Define an HyperCube/Simplex class
setClass("HCSimplex", representation(
	b = "array",
	h = "numeric",
	s = "numeric",
	iteration = "numeric"
	)
)

setMethod("initialize", "HCSimplex", function(.Object, or, lim, multiplier, it) {
  .Object@b <- or
	.Object@h <- lim
	.Object@s <- multiplier
	.Object@iteration <- it
	.Object
})

VHmeasure <-function(frame, hc) {
	if (class(frame) != "AEP.frame") 
		stop("The parameter 'frame' is not of class 'AEP.frame'")
	if (length(hc@b) != frame@nb_var) {
		stop("hc@b is too long")		
	}
	samplePoints <- array(0, c(frame@N, frame@nb_var))
	Uval <- array(0, c(frame@N, frame@nb_var))
	for (i in 1:frame@N) {
		samplePoints[i,] <- hc@b + (hc@h * frame@binary_matrix[i,])
		for (j in 1:frame@nb_var) Uval[i,j] <- frame@ECDF(j, samplePoints[i,j])	
	}
	H <- C.n(u=Uval, U=frame@pobs_data)
	Hm <- frame@Hmultiplier * H	
	sum(Hm)
}

AEP <- function(frame, data, value, alpha, nb_iteration) {
  if (missing(frame)) {
    if (missing(data) || missing(alpha) || missing(nb_iteration)) stop(AEP.error.input_msg) 
    if (class(data) != "data.frame") stop(AEP.error.input_msg)    
    frame <- new("AEP.frame", data, alpha, nb_iteration)    
  } else {
    if (!(missing(data) & missing(alpha) & missing(nb_iteration))) stop(AEP.error.input_msg)     
    if(class(frame) != "AEP.frame") stop(AEP.error.input_msg)
  }
  resolveAEP(frame, value)
}

resolveAEP <- function(frame, value) {
#	nb_of_resolution <- 0
	result <- 0
	spx <- new("HCSimplex", array(0, c(frame@nb_var)), value, 1, 1)
	ListOfSimplexes <- c(spx)	
	while (length(ListOfSimplexes) > 0) {
		spx <- ListOfSimplexes[[1]]
		if (length(ListOfSimplexes) > 1) ListOfSimplexes <- ListOfSimplexes[2:length(ListOfSimplexes)]
		else ListOfSimplexes <- NULL
		hc <- new("HCSimplex", spx@b, frame@alpha * spx@h, spx@s, spx@iteration)
		if (spx@iteration < frame@nb_iteration) {
			new.it <- spx@iteration + 1
			for (j in 2:frame@N) {
				new.b <- spx@b + (frame@alpha * spx@h * frame@binary_matrix[j,])
				new.h <- (1 - (frame@alpha * frame@sum_ones[j])) * spx@h
				new.s <- spx@s * frame@m[j]				
				ListOfSimplexes <- c(ListOfSimplexes, new("HCSimplex", new.b, new.h, new.s, new.it))			
			}
		}
		result <- result + hc@s * VHmeasure(frame, hc)
#		nb_of_resolution <- nb_of_resolution + 1
	}
#	print(nb_of_resolution)
	result	
} 

       
			
test.debug <- FALSE

if (test.debug) {
  
	#test case
	nb.obs <- 30000
	nb.var <- 2
	mydata = matrix(runif(nb.obs * nb.var), c(nb.obs, nb.var))
	mydata <- mydata * 50
  mydata <- data.frame(mydata)
  
	my.aep.frame.test <- new("AEP.frame", data=mydata, alpha=0.75, nb_iteration=4)
	print(AEP(frame=my.aep.frame.test, value=25))
	print(AEP(frame=my.aep.frame.test, value=50))
  
  print(AEP(data=mydata, value=25, alpha=0.75, nb_iteration=4))
	print(AEP(data=mydata, value=50, alpha=0.75, nb_iteration=4))
  
	#end of test case
  
} else {	
  
  par(ask=TRUE)
	par(mfrow=c(1,1))
	mydata <- read.table("output_10_100.txt", header=TRUE, sep=" ")
	observed.exectime <- array(0, c(dim(mydata)[1]))
	for (i in 1:length(observed.exectime)) observed.exectime[i] <- sum(mydata[i,])
	plot( x = 1:length(observed.exectime), y = observed.exectime, 
		main = "Observed execution time",
		xlab = "#run",
		ylab = "Observed Exec Time",
		type = "l", 
		col = "red")

	observed.ecdf <- ecdf(observed.exectime)
	plot(observed.ecdf,
		main = "Observed ECDF",
		xlab = "#run",
		ylab = "P[ExecTime < X]"
	)
  
  	
	my.step <- (max(observed.exectime) - min(observed.exectime)) / 100
	estimated.ecdf <- array(0, c(100))

	my.aep.frame <- new("AEP.frame", data=mydata, alpha=0.5, nb_iteration=2)
  
	index <- 1
	steps <- seq(min(observed.exectime), 2*max(observed.exectime), by=my.step)
	for (i in steps) {
		estimated.ecdf[index] <- AEP(frame=my.aep.frame, value=i)
		index <- index + 1
	}
	
	plot( x = steps, y = estimated.ecdf, 
		main = "Estimated CDF",
		xlab = "Estimated CDF",
		ylab = "P[ExecTime < X]",
		type = "l", 
		col = "red")
}
